This notebook is part of PyBroMo a python-based single-molecule Brownian motion diffusion simulator that simulates confocal smFRET experiments.
In [ ]:
%matplotlib inline
import numpy as np
import tables
import matplotlib.pyplot as plt
import seaborn as sns
import pybromo as pbm
print('Numpy version:', np.__version__)
print('PyTables version:', tables.__version__)
print('PyBroMo version:', pbm.__version__)
As a memo, let's write some formulas related to the FRET efficiency:
$$ k = \frac{F_a}{F_d} \quad,\qquad E = \frac{k}{k+1} \qquad\Rightarrow\qquad k = \frac{E}{1-E}$$
In [ ]:
S = pbm.ParticlesSimulation.from_datafile('016', mode='w')
In [ ]:
def em_rates_from_E(em_rate_tot, E_values):
E_values = np.asarray(E_values)
em_rates_a = E_values * em_rate_tot
em_rates_d = em_rate_tot - em_rates_a
k_values = E_values/(1 - E_values)
assert np.allclose((em_rates_a/em_rates_d), k_values)
em_rates = np.hstack([em_rates_a, em_rates_d])
return em_rates
In [ ]:
em_rate_tot = 200e3
E_list = np.array([0.01, 0.02, 0.05, 0.1, 0.2, 0.4])
em_rate_list = em_rates_from_E(em_rate_tot, E_list)
em_rate_list
In [ ]:
# Get the random state at the end of the diffusion simulation
saved_rs_state = S.traj_group._v_attrs['last_random_state']
pbm.hash_(saved_rs_state)
In [ ]:
em_rate_list
Simulate timestamps for background = 1kcps:
In [ ]:
rs = np.random.RandomState()
rs.set_state(saved_rs_state)
In [ ]:
%%timeit -n1 -r1
for em_rate in em_rate_list:
print(' Emission rate: ', em_rate, flush=True)
S.simulate_timestamps_mix(max_rates=(em_rate,), populations=(slice(0, 20),),
bg_rate=1e3, rs=rs)
Simulate timestamps for background = 4kcps:
In [ ]:
%%timeit -n1 -r1
for em_rate in em_rate_list:
print(' Emission rate: ', em_rate, flush=True)
S.simulate_timestamps_mix(max_rates=(em_rate,), populations=(slice(0, 20),),
bg_rate=4e3, rs=rs)
In [ ]:
for k in S.ts_store.h5file.root.timestamps._v_children.keys():
if not k.endswith('_par'):
print(k)
In [ ]:
ts, ts_par = S.get_timestamps_part('Pop1_P20_Pstart0_max_rate198000cps_BG4000cps_t_1s_rs_8798a6')
In [ ]:
ts[:]
In [ ]:
bins = np.arange(0, 1, 1e-3)
plt.hist(ts*ts.attrs['clk_p'], bins=bins, histtype='step');
In [ ]:
group = '/timestamps'
print('Nodes in in %s:\n' % group)
print(S.ts_store.h5file.get_node(group))
for node in S.ts_store.h5file.get_node(group)._f_list_nodes():
print('\t%s' % node.name)
#print('\t %s' % node.title)
In [ ]:
[t for t in S.timestamp_names if 'BG4000cps' in t]
In [ ]:
S.ts_store.close()
In [ ]:
In [ ]: